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Abstract — This  paper  presents  a selective  wavelet  shrinkage 
algorithm  for  digital  image  denoising.  The  performance  of  this 
method  is  an  improvement  upon  other  methods  proposed  in 
the  literature  and  is  algorithmically  simple  for  large  computa- 
tional savings.  The  improved  performance  and  computational 
speed  of  the  proposed  wavelet  shrinkage  algorithm  is  pre- 
sented and  experimentally  compared  with  established  methods. 
The  denoising  methodology  incorporated  in  this  new  algorithm 
involves  a two-threshold  validation  process  for  real-time  selec- 
tion of  wavelet  coefficients.  The  two-threshold  criteria  selects 
wavelet  coefficients  based  on  their  absolute  value,  spatial  reg- 
ularity, and  regularity  across  multiresolutional  scales. 

The  proposed  algorithm  takes  image  features  into  considera- 
tion in  the  selection  process.  Statistically,  most  images  have 
regular  features  resulting  in  connected  subband  coefficients. 
Therefore,  the  resulting  subbands  of  wavelet  transformed  im- 
ages in  large  part  do  not  contain  isolated  coefficients.  In  the 
proposed  algorithm,  after  coefficients  are  selected  due  to  their 
magnitude,  image  features  in  terms  of  spatial  regularity  are 
used  to  further  reduce  the  number  of  coefficients  kept  for  im- 
age reconstruction. 

The  proposed  wavelet  denoising  technique  is  unique  in  that  its 
performance  improved  upon  several  other  established  wavelet 
denoising  techniques  as  well  as  being  computationally  efficient 
to  facilitate  realtime  image  processing  applications. 

1.  Introduction 

The  recent  advancement  in  multimedia  technology  has  pro- 
moted an  enormous  amount  of  research  in  the  area  of  im- 
age and  video  processing.  Included  in  the  many  image  and 
video  processing  applications  such  as  compression,  enhance- 
ment, and  target  recognition  are  preprocessing  functions  for 
noise  removal.  Noise  removal  is  one  of  the  most  common  and 
important  processing  steps  in  many  image  and  video  systems. 

Because  of  the  importance  and  commonality  of  preprocessing 
in  most  image  and  video  systems,  there  has  been  an  enormous 
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amount  of  research  dedicated  to  the  subject  of  noise  removal, 
and  many  different  mathematical  tools  have  been  proposed. 
Variable  coefficient  linear  filters  [8, 10],  adaptive  median  fil- 
ters [3],  DCT  based  solutions  [9],  cluster  filtering  [11],  etc. 
have  all  been  proposed  in  the  literature. 

The  wavelet  transform  has  also  been  used  to  suppress  noise 
in  digital  images.  It  has  been  shown  the  reduction  of  abso- 
lute value  in  wavelet  coefficients  is  successful  in  signal  restora- 
tion [5].  This  process  is  known  as  wavelet  shrinkage.  Other, 
more  complex  denoising  techniques  select  or  reject  wavelet 
coefficients  based  on  their  predicted  contribution  to  recon- 
structed image  quality  . This  process  is  known  as  selective 
wavelet  shrinkage,  and  many  works  have  used  it  as  the  pre- 
ferred method  of  image  denoising  [2,4-7]. 

Two  of  the  frontrunners  in  selective  wavelet  shrinkage  for  the 
removal  of  noise,  Mallat  and  Hwang  prove  the  successful  re- 
moval of  noise  in  signals  via  the  wavelet  transform,  by  select- 
ing and  rejecting  wavelet  coefficients  based  on  their  Lipschitz 
exponents  [5].  Although  this  fundamental  work  in  image  de- 
noising is  successful  in  the  removal  of  noise,  its  application  is 
broad  and  not  focused  on  image  noise  removal,  so  the  results 
are  not  optimal. 

Malfait  and  Roose  refined  the  selective  shrinkage  denoising 
approach  by  applying  a Bayesian  probabilistic  formulation, 
and  modelled  the  wavelet  coefficients  as  Marcov  random  se- 
quences [6].  This  method  is  focused  on  Image  de-noising  and 
its  results  are  an  improvement  upon  [5]. 

Later,  Pizurica,  et  al.  continued  on  the  work  done  by  [6].  This 
work  applied  a statistical  probabilistic  model  that  extended  to 
both  spatial  and  multiresolutional  coefficient  redundancies  to 
decide  on  important  and  unimportant  values  [7]. 

Although  both  algorithms  in  [6]  and  [7]  give  adequate  results 
in  denoised  image  quality,  their  computational  complexities 
make  them  impractical  for  most  image  and  video  processing 
applications. 

We  have  developed  a new  selective  wavelet  shrinkage  algo- 
rithm which  has  a distinct  aspect.  It  uses  a double  validation 
process  to  select  wavelet  coefficients  for  quality  image  recon- 
struction. Furthermore,  the  algorithm  is  much  more  efficient  in 
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computational  complexity  than  previous  works.  The  algorithm 
is  comprised  of  three  processing  steps  which  are  common  to  all 
selective  wavelet  shrinkage  algorithms  mentioned  previously. 
First,  a corrupted  image  is  decomposed  into  multiresolutional 
subbands  via  the  wavelet  transform.  Next,  wavelet  coefficients 
are  selected  or  rejected  based  upon  criteria  developed  by  the  al- 
gorithm designer.  Finally,  the  denoised  image  is  formed  by  re- 
constructing the  remaining  coefficients  via  the  inverse  wavelet 
transform.  The  processing  step  of  most  cost  computationally 
and  most  important  in  denoising  performance  is  the  coefficient 
selection  process,  which  calls  for  effective  and  efficient  criteria 
to  select  or  reject  wavelet  coefficients. 

This  double  validation  method  for  selecting  wavelet  coeffi- 
cients results  in  a denoising  algorithm  which  gives  improved 
results  upon  [5-7],  but  without  the  computational  complexity. 
The  two-threshold  requirement  investigates  the  regularities  of 
wavelet  coefficients  both  spatially  and  across  scales  for  predic- 
tive coefficient  selection,  providing  selective  wavelet  shrinkage 
to  non-decimated  wavelet  subbands. 

Following  the  Introduction,  Section  2 gives  a wavelet  overview 
and  theory  on  the  2D  non-decimated  wavelet  analysis  and 
synthesis  filters.  The  proposed  denoising  algorithm  is  ap- 
plied to  non-decimated  wavelet  coefficients.  Section  3 then 
describes  the  coefficient  selection  process  prior  to  selective 
wavelet  shrinkage.  Section  4 gives  testing  results  for  param- 
eter selection.  Section  5 gives  the  estimation  algorithms  for 
proper  parameter  selection,  and  Section  6 gives  the  results  and 
discussion. 

2.  Wavelet  Overview 

Wavelet  Filterbank  Theory 

Let  ak  [n]  be  scaling  coefficients  of  scale  k and  position  n,  and 
let  h[n ] be  the  filter  coefficients  corresponding  to  the  scaling 
function.  From  wavelet  theory,  we  know 

ak+1  [n]  = '^2ak[m]h[m  - 2n],  (1) 

771 

where  the  coefficients  ak+l[n]  represent  a coarser  resolution 
then  ak[n).  Equation  1 indicates  that  the  scaling  function  co- 
efficients ak+1  scale  can  be  obtained  by  convolving  a reversed 
h\n\  with  ak,  and  downsampling  by  two. 

Very  similarly, 

4+1  M = ak[m}g[m  - 2n],  (2) 

m 

where  dk[n]  are  the  wavelet  coefficients  of  scale  k and  position 
n.  g[n]  is  the  set  of  filter  coefficients  corresponding  to  the 
wavelet. 

From  Equations  I and  2,  we  can  obtain  increasingly  coarser 
scales  of  scaling  coefficients,  ak+1,  and  wavelet  coefficients. 


dk+  j [n],  by  convolving  the  scaling  function  coefficients,  o,k  [n] 
by  both  a reversed  scaling  function  filter,  h[n],  a reversed 
wavelet  filter,  g[n],  and  downsampling  by  two.  Figure  1 gives 
a block  diagram  of  the  wavelet  analysis  filterbank. 


Figure  1 - Wavelet  Decomposition 


Because  each  filtered  output  is  downsampled  by  two,  the  same 
number  of  total  coefficients  remains  the  same,  regardless  of  the 
number  of  resolution  levels,  k. 

The  reconstruction  of  finer  scaling  coefficients  is  obtained  by, 

«fc[n]  = Em  ak+ 1 M h[n  - 2m] 

+ i2mdk+1[™]g[n  - 2m]. 

From  Equation  3,  we  can  the  obtain  an  arbitrarily  fine  scale 
representation  of  a signal  by  upsampling  the  scaling  and 
wavelet  coefficients,  and  filtering  the  coefficients  with  their 
respective  filters,  h[n]  and  g[n].  The  wavelet  reconstruction 
block  diagram  is  given  in  Figure  2. 


Figure  2 - Wavelet  Reconstruction 


Non-Decimated  Wavelet  Transform 

In  certain  applications  such  as  signal  denoising,  it  is  not  de- 
sirable to  downsample  the  wavelet  coefficients  after  decom- 
position, because  the  spatial  resolution  of  the  coefficients 
is  degraded  due  to  downsampling.  Therefore,  for  the  non- 
decimated  case,  each  subband  contains  the  same  number  of 
coefficients  as  the  original  signal.  So  we  must  have 

c*fc[2fc+1n]  = ak[n\ 

\k[2k+1n]  = dk[n], 

where  otk[n]  are  the  non-decimated  scaling  function  coeffi- 
cients, and  Afc  [n]  are  the  non-decimated  wavelet  coefficients. 
We  can  substitute  Equation  4 into  Equation  1 to  find 

ak+ 1 N = Em  h[m]ak[m  - 2n] 

ak+i  [2 fc+2n]  = Em  h[m]ak[2k+1(m  - 2 n)]  (5) 

afe+iEn]  = Em  h[m\ak[2k+1m  - n]. 
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The  2fc+1  scalar  introduced  into  Equation  5 is  the  equivalent  of 
upsampling  h[n ] by  k + 1 prior  to  its  convolution  with  [n] . 
Similarly  we  substitute  Equation  4 into  Equation  2 to  obtain, 

Afc+1  [n]  = ^raff[mK[2H1m  - n].  (6) 

Figure  3 gives  a block  diagram  of  the  non-decimated  wavelet 
decomposition. 


where 

f(x,y)  = otu-i[x,y],  (13) 

The  four  coefficent  sets  given  in  Equation  12  is  referred  to  as 
the  low-low  band,  an <k+l , the  high-low  band,  A Wjfc+1 , the  low- 
high  band,  Xih,k+1 , and  the  high-high  band,  A hh,k+t  ■ The  sub- 
bands are  named  due  to  the  order  in  which  the  scaling  and/or 
the  wavelet  filters  process  the  scaling  function  coefficients. 


Figure  3 - Non-decimated  wavelet  decomposition 


The  synthesis  of  the  non-decimated  wavelet  transform  also  dif- 
fers from  the  downsampled  case.  From  Equation  3,  we  have 

ak[2n]  = ]T/i[2(n  - m)]afc+1  [m]  + s[2(n  - m)\dk+l  [m], 

m m 

(7) 

Substituting  (p  = n - m)  we  obtain, 
ak[2n]  = ^/i[2p]afc+1[n-p]  + ^ff[2p]4+1[n-p],  (8) 

P V 

Substituting  Equation  4 into  Equation  8 we  have, 

<*fc[2fc+2n]  = EP  h[2p]ak+l  [2 k+2(n  - p)] 

+ p[2p]Afc+1  [2fc+2(n  — p)]  ’ 

and 

akin)  = EP  h{2p\ak+l  [n  - 2fc+2p] 


(9) 


+ Ept?[2p]Afc+1[n-2fc+2p]. 


(10) 


Looking  at  Equation  10  information  is  being  thrown  away  by 
downsampling  ak+J  [n]  and  Xk+i  [n]  by  2 prior  to  convolution. 
Because  the  downsampling  in  the  analysis  filters  is  eliminated, 
a downsample  by  2 is  shown  in  the  synthesis  equation,  Equa- 
tion 10.  If  a downsample  by  2 is  not  performed,  i.e.  (m  = 2p), 
then  we  must  divide  by  2 to  provide  power  equality.  That  is. 


ak[n } = \ Em  hlm}ak+!  ln  - 2fc+1m] 
+]  Em5HAfc+1[n  - 2fc+1m] 


(11) 


Figure  4 gives  a block  diagram  of  the  non-decimated  wavelet 
transform  synthesis. 


KM 


«»*>! 


For  the  synthesis  of  / we  have 

| Em,n  ^[m]^[n]a«,*;+1  lx  ~ 2fc+1m,  y — 2fe+1n] 
+ 5 E m,n  h[m}g[n\\hl'k+1  [x  - 2 k+1m,  y - 2 fc+1n] 
+ 1 Em,n  9[m]h[n}\lh'k+1  [x  - 2 k+1m,  y - 2k+1n j 
+ 3 Em  - 2fc+1m,p  - 2 fc+1n] 

(14) 

3.  Coefficient  Support 

One  of  the  many  advantages  of  the  wavelet  transform  over 
other  mathematical  transformations  is  the  retention  of  spatial 
information  in  the  wavelet  domain.  Because  of  this  informa- 
tion retention,  there  exists  a spacial  regularity  in  the  sub-bands 
of  wavelet  transformed  images.  Statistically,  most  images  have 
regular  features  resulting  in  connected  subband  coefficients. 
Therefore,  the  resulting  subbands  of  wavelet  transformed  im- 
ages in  large  part  do  not  contain  isolated  coefficients.  This 
regularity  can  aid  in  deciding  which  coefficients  should  be  se- 
lected for  reconstruction,  and  which  should  be  thrown  away 
for  maximum  reconstructed  image  quality.  The  correlation  be- 
tween coefficients  in  wavelet  sub-bands  has  been  discovered 
by  many  works,  but  our  method  in  which  this  spacial  regular- 
ity is  exploited  is  unique. 

We  start  with  an  image  signal  corrupted  with  additive  noise, 
i.e. 

f{x,y)  = f{x,y)  + r](x,y),  (15) 

where  f(x,y ) is  the  noiseless  2D  signal,  rj(x,y)  is  a random 
noise  function,  and  f(x,  y)  is  the  corrupted  signal. 

The  wavelet  transform  of  f(x,y)  generates  coefficients,  X.ik 
using  Equations  12  and  13.  X.tk  is  used  to  create  a boolean 
coefficient  map,  I.  k. 

0® 


Figure  4 - Non-decimated  wavelet  synthesis 

We  can  expand  to  the  two-dimensional  case.  For  a 2-D  discrete 
signal,  / we  have, 

au<k+l[x,y}  = J2n,m  h{n}h{Tn}aiitk[2k+1ni  - x,2k+1n  - y] 

A w,fc+1  [ar, y\  = ]Cn,m  h{n}g[m}aink[2k+1m  - x,2 k+1n-  y] 
Ai/a+ikyj  = Y.n,m9ln]hlm]aiiM{2k+lm  - x,2 k+1n  - y) 
Ah/,^,^,^]  = Y2nm9ln]glrn]autk{2k+1m-x,2k+ln-y), 

(12) 


A valid  coefficient  is  defined  as  a coefficient  value,  A.^x,  y], 
which  resuls  in  I.tk[x,y\  = 1,  hence  the  coefficient  has  been 
selected  due  to  its  magnitude. 

After  coefficients  are  selected  by  magnitude,  spacial  regularity 
is  used  to  further  reduce  the  number  of  coefficients  kept  for 
image  reconstruction. 

From  I.>k  we  can  count  the  number  of  support  pixels  around 
a particular  I.ik[x,y],  S.ik\x, y]  is  the  sum  of  all  I.)k  which 
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support  the  current  boolean  value  I.tk[x,y];  that  is,  the  total 
number  of  all  valid  coefficient  values  which  are  spatially  con- 
nected to  y\. 


Not  only  do  coefficients  of  wavelet  transformed  images  have 
spatial  regularity,  they  also  have  regularity  across  scales. 
Therefore,  to  exploit  this  regularity,  we  have 


We  say  that  a coefficient  is  spatially  connected  to  another  if 
there  exists  a path  of  valid  coefficients  between  the  two,  in 
any  direction.  Figure  5 gives  a generic  coefficient  map.  The 
valid  coefficients  are  highlighted  in  gray.  From  Figure  5 it  can 
be  shown  that  coefficients  A,  B,  C,  and  H do  not  support  any 
other  coefficients  in  the  coefficient  map.  However,  coefficients 
D and  F support  each  other,  coefficients  E and  G support  each 
other,  and  N and  O support  each  other.  Also,  coefficients  I,  J, 
K,  L,  M,  P,  and  Q all  support  one  another. 


Figure  S - Generic  Coefficient  Array 


Figure  6 gives  the  value  of  S.tk  [x,  y]  for  each  of  the  valid  coef- 
ficients given  in  Figure  5. 


Figure  6 - Generic  Coefficient  Array,  with  corresponding  S.  *.  values 


An  algorithm  for  computing  S.^rr,  y]  is  given  in  the  Ap- 
pendix. S.tk  is  used  to  refine  the  original  coefficient  map  J.  fc 
by 


J-,k[x^y\  = { 


l.&fcfc.J/]  > s 
0,  else 


(17) 


where  J.,fc[x,  y)  is  the  refined  coefficient  map,  and  s is  the  nec- 
essary number  of  support  coefficients. 


L,k 


[x,  y]  = { 


0, else 

(18) 


The  de-noised  image  is  then  reconstructed  using  the  supported 
coefficient  values,  L fc[x,y]  in  the  synthesis  equation  given  in 
Equation  14.  Thus,  we  have 

y]  = 5 T,m,n  h[m\h[n}aii,k+1  [z  - 2 k+1m,  y - 2fc+1n] 

+ \ Em,n  h[m]g\n]Lhi)k+1  [x  - 2fc+1m,  y - 2fc+1n] 

+ 3 Em, n 9[m]h[n]Lih,k+1  [x  - 2 k+1m,  y - 2fc+1n]  ’ 
+ 3 Em,n  9[m\g [n]Lhh,k+1  [x  - 2*+1m,  y - 2fc+1n] 

(19) 

where 

f(x,y)  = au-  \[x,y].  (20) 

dn.fc  and  / are  the  reconstructed  scaling  function  coefficients 
and  denoised  image,  respectively. 


4.  Test  Images  and  Selection  of 

THRESHOLD  r AND  SUPPORT  s 

In  the  above  algorithm,  we  select  wavelet  coefficient  values 
based  upon  a threshold  value  r and  a support  value  .s.  So  now 
we  must  obtain  choices  for  these  values  for  optimal  image  de- 
noising.  We  start  with  a series  of  test  images,  given  in  Figure 
7.  These  test  images  are  to  be  used  in  choosing  the  values  for 
r and  s. 


Figure  7 - Test  Images 


The  test  images  are  all  256x256  in  size.  Starting  from  the  up- 
perleft  image  and  going  clockwise,  the  images  are  ’’Lenna”, 
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’’Airplane”,  ’’Fruits”,  and  ’’Girl”.  All  of  these  images  are  well 
known  standards  in  the  area  of  image  processing. 


We  begin  testing  for  the  optimal  values  of  r and  s by  artificially 
adding  Gaussian  noise  to  each  of  the  four  images,  denoising  all 
four  images  with  a particular  r and  s,  and  recording  the  aver- 
age Peak  Signal-to-Noise  Ratio  (PSNR).  PSNR  of  an  image  is 
defined  by 


PSNR  = 20  log10 


255  \ 
y/mse])  ’ 


(21) 


where 

msef  = WN  ^ ^ ~ Kx'  • (22) 


PSNR  Results  with  noise,  o s 30 


spatial  support  pixels,  s 


Threshold  value,  t 


M and  N are  the  width  and  height  of  the  images,  respectively. 


Figure  8 - PSNR  Results  for  Test  images,  an  = 30 


The  Haar  wavelet  is  selected  for  image  denoising: 


, n = 0, 1 
,else 


g[n]  = { 


-1 


y/2 

0 


,n  = 0 
,n  = 1 
, else 


(23) 


The  Haar  wavelet  is  used  in  a non-decimated  wavelet  decom- 
position of  the  original  image.  6 subband  levels  are  used,  i.e. 
k = —1  to  5.  The  proposed  selective  wavelet  shrinkage  algo- 
rithms is  applied  to  all  wavelet  subbands,  and  the  subbands  are 
synthesized  by  the  non-decimated  wavelet  inverse  transform. 


Preliminary  tests  had  shown  the  Haar  wavelet  has  the  most 
promise  in  reconstructed  image  quality.  The  compact  support 
of  the  Haar  wavelet  enables  the  generation  of  wavelet  coeffi- 
cients which  represent  the  least  amount  of  original  pixel  data. 
Therefore,  when  a coefficient  is  removed  because  of  its  in- 
significance. The  result  affects  the  least  number  of  pixels  in 
the  reconstruction. 


We  recorded  each  of  the  PSNR  averages  for  r ranging  from 
0 — 100  and  s ranging  from  0 — 15.  We  tested  the  proposed 
algorithm  by  applying  additive  white  Gaussian  noise  (AWGN) 
with  a standard  deviation  (cr„)  of  10, 20, 30, 40,  and  50,  to  each 
of  the  test  images.  Our  method  of  selective  wavelet  shrinkage 
is  applied  to  the  corrupted  image,  and  the  resulting  PSNR  is 
recorded.  The  results  of  the  testing  in  which  an  = 30  is  given 
in  Figure  8. 

Figure  9 gives  the  r and  s which  provide  the  largest  average 
PSNR  for  each  noise  level.  We  will  refer  to  these  particular 
values  as  rm  and  sm.  Figure  9 suggests  that  parameters  rm 
and  sm  are  functions  of  the  standard  deviation  of  the  artificial 
noise,  an. 


5.  Estimation  of  Parameter  Values 

Noise  Estimation 


Noise  Level  (crn) 

10: 

20: 

30 

40 

50 

Max,  Avg  PSNR 

33.42 

30.10 

28  34 

27  07 

26.14 

Sm  value 

2 

3 

5 

9 

11 

Hn  value 

14 

28 

38 

46 ; 

62 

Figure  9 - Maximum  average  PSNR  of  test  images  for  various  noise  levels 
and  their  corresponding  threshold  and  support  values 


tain  an  estimate  of  the  optimal  values  for  r and  s from  the 
standard  deviation  of  the  noise  level.  However,  the  level  of 
noise  in  a given  digital  image  is  unknown.  So  we  must  first 
estimate  the  noise.  Several  well  known  algorithms  have  been 
given  in  the  literature  to  estimate  image  noise.  From  [1,7]  a 
median  value  of  the  A hh,o  subband  is  used  in  the  estimation 
process.  However,  retrieving  the  median  value  in  a subband 
requires  a sorting  algorithm  which  is  computationally  expen- 
sive. We  propose  an  averaging  noise  estimation  algorithm.  A 
threshold  value,  e is  used  to  sort  out  strong  signal  (edge)  coef- 
ficients from  the  noisy  data. 

I-WMII.  (24) 

x,y 

where  e is  a measure  of  the  average  magnitude  of  a coefficient 
in  the  Xhh,o  subband.  We  also  have. 


(25) 


In  our  noise  estimation,  we  use  an  average  of  wavelet  coef- 
ficients instead  of  a median.  The  result  is  a computationally 
simpler  estimate  of  the  image  noise.  p[. r,  y\  is  used  as  a refine- 
ment parameter  to  remove  large  signal  values  located  in  the 
Xhh,o  subband.  We  estimate  the  noise  by, 

2 r 

= -^^2p[x,y}\Xhhfi[x,y]\  (26) 


It  can  be  shown  from  the  values  given  in  Figure  9 that  the  pa- 
rameters rm  and  sm  are  functions  of  a„,  therefore  we  can  ob- 


where  K is  the  number  non-zero  terms  in  the  summation  of 
Equation  26. 
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Parameter  Estimation 


6.  Results 


Using  the  known  level  of  noise  added  to  the  original  images, 
we  can  estimate  the  values  of  Tm  and  sm  given  in  Figure  9. 
We  use  the  LMMSE  (Linear  Minimum  Mean  Squared  Error) 
as  the  estimation  procedure.  That  is,  we  find  two  parameters 
aT  and  bT  such  that 

Tmi^n)  — ar(Tn  + bT.  (27) 

The  choise  aT  and  bT  will  minimize  the  mean  squared  error 
(mse): 

mseT  = ^(rm(<r„)  - f^(crn))2.  (28) 

On 

Similarly,  we  find  an  estimate  of  s 

7 -f-  bs  (29) 

where  as  and  bs  are  chosen  to  minimize 

mses  = ^(sm(u„)  - S^K))2.  (30) 

On 


From  the  data  given  in  Figure  9,  the  values  of  aT, 
bs.  are  found 

aT  = 1.14 
bT  = 3.40 
as  = 0.21 
bs  = 0.60 


Qj,  and 


(31) 


The  values  of  rm  and  sm  are  given  in  Figure  10  as  well  as  the 
corresponding  LMMSE  estimates,  given  in  Equation  31.  As 
given  in  Figure  10  the  estimated  values  are  the  best  linear  fit 
into  the  data.  Note  that  the  support  value  s must  be  an  integer 
value. 


threshold  and  support  estimation  based  upon  noise  level 
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Figure  10  - Optimal  values  for  rm  and  sm,  and  their  corresponding  esti- 
mates, TVn  and  s^i 


The  images  ’’peppers”  and  ’’house”  are  used  for  gauging  the 
performance  of  our  denoising  algorithm.  These  two  images 
have  also  been  used  in  the  results  of  [5-7].  Therefore,  we  can 
compare  our  performace  with  other  recent  algorithms  given  in 
the  literature. 

We  have  corrupted  both  the  peppers  image  and  house  image 
AWGN  and  used  the  proposed  method  for  denoising.  The  re- 
sults are  given  in  Figures  12  and  13. 
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Figure  11  - PSNR  Comparison  of  the  of  the  proposed  method  to  other  meth- 
ods in  the  literature 


Figure  12  - Results  of  the  proposed  image  denoising  algorithm.  Top  left: 
Original  ’’peppers”  image,  top  right:  Corrupted  Image,  c„  = 37.75  - PSNR 
= 16.60,  bottom  left:  denoised  image  with  estimated  r and  s = 0 - PSNR  = 
25.76  , bottom  right:  denoised  image  with  estimated  r and  s - PSNR  = 27.05 


The  threshold  value  t and  the  support  value  s will  be  deter- 
mined by 

T = aTOn  + bT 

s = asan  + bs 


Figure  1 1 gives  the  results  of  our  proposed  method,  as  well 
as  the  results  of  [5-7].  As  shown  in  Figure  11,  the  results  of 
the  proposed  method  are  an  improvement  to  other  methods  de- 
scribed in  the  literature. 
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Figure  13  - Results  of  the  proposed  image  denoising  algorithm.  Top  left: 
Original  ’’house”  image,  top  right:  Corrupted  Image,  ern  = 32.47  - PSNR  = 
17.90,  bottom  left:  denoised  image  with  estimated  r and  s = 0 - PSNR  = 
28.10  , bottom  right:  denoised  image  with  estimated  r and  s - PSNR  = 30.08 

In  addition  to  improved  performance,  the  proposed  algorithm 
is  computationally  simple  to  facilitate  real-world  applications. 
However,  computational  simplicity  is  a fairly  difficult  metric 
to  determine  when  dealing  with  older  literature.  The  process- 
ing ability  of  modem  computers  gives  more  recent  literatures’ 
algorithms  processing  time  an  unfair  advantage.  However,  the 
computation  time  of  the  proposed  method  is  at  least  an  order 
of  magnitude  greater  than  previous  methods,  and  we  have  done 
testing  on  some  older  machines  for  a more  accurate  compari- 
son. Figure  14  gives  the  computational  results  of  the  proposed 
method  as  well  as  the  results  of  [6,7]. 


: Processor 
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^Proposed  Alqorithm 
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I Pizurica  3-band.  [6] 

45.00 

*** 

iPizurica  2-band,  [6| 

30.00 
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iMalfait  and  Roose,  [4] 

**-*•  **-» 

180.00 

Figure  14  - Computation  Times  for  a 256x256  Image,  in  seconds 


Although  this  is  not  a true  comparison  of  the  difference  com- 
putational complexity  between  the  proposed  algorithm  and  that 
of  [6,7],  the  proposed  algorithm  does  show  a substantial  drop 
in  computation  time. 

In  this  paper,  a new  selective  wavelet  shrinkage  algorithm  for 
image  denoising  has  been  described.  The  proposed  algorithm 
uses  a two-threshold  support  criteria  which  investigates  coef- 
ficient magnitude,  spatial  support,  and  support  across  scales 
in  the  coefficient  selection  process.  The  computationally  sim- 
ple algorithm  facilitates  real  world  applications  and  the  perfor- 
mance results  are  an  improvement  upon  established  methods 


described  in  the  literature. 


7.  Appendix 

The  computation  of  S.^[x,y\  if  given  from  the  following  al- 
gorithm: 

N{)  = {[-1,-1],  [-1,0],  [-1,1],  [0,-1], 

[0,1],  [1,-1],  [1,0],  [1,1]} 

0 = 0,  t = 0,  p = 0,  l).,fc(0)  = (x, y) 
if  I-,k[x,y]=—  1, 

while  D.tk{t)  7^  NULL, 

( i,j ) = D;k(t) 
t = t + 1 

for  m = 0 to  7, 

if  ((I,k[(iJ)+N(m)]==  1)  (33) 

and  ( 0[(i,j ) + N(m)\  ==  0)), 

p = p+  1 

D,k(p)  = ((iJ)  + N(m)) 

0[(i,j)  + N(m)}  = 1, 
end  if 
end  for 
end  while 
end  if 

S.,k[x,y]  = t 

0[x,y]  is  a boolean  value  to  determine  whether  a particular 
J.,fc[x,  y]  value  has  been  counted  previously.  D is  an  array  of 
spatial  coordinates  of  valid  coefficients  that  support  the  current 
coefficient  value  I.tk  [x,  y].  N is  a set  of  vectors  corresponding 
to  neighboring  coefficient  values. 
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